library(Matrix)
library(aroma.light)
## Registered S3 method overwritten by 'R.oo':
## method from
## throw.default R.methodsS3
## aroma.light v3.14.0 (2018-09-04) successfully loaded. See ?aroma.light for help.
countsFull <- readMM("~/Downloads/GSE128889_RAW/GSM3717977_SCmurinep12_matrix.mtx.gz")
ft <- read.table("~/Downloads/GSE128889_RAW/GSM3717977_SCmurinep12_genes.tsv.gz")
bc <- read.table("~/Downloads/GSE128889_RAW/GSM3717977_SCmurinep12_barcodes.tsv.gz")
rownames(countsFull) <- as.character(ft[,2])
keep <- rowSums(countsFull > 1) >= 400
counts <- countsFull[keep,]
rownames(counts) <- ft[keep,2]
countsFQ <- normalizeQuantileRank(as.matrix(counts))
rownames(countsFQ) <- ft[keep,2]
pca <- prcomp(t(log1p(countsFQ)), scale. = FALSE)
rd <- pca$x[,1:8]
plot(rd[,1:2], pch=16, cex=1/2)
plotByGene <- function(rd, geneCount, ng=10, main=NULL, ...){
pal <- wesanderson::wes_palette("Zissou1", n=ng, type="continuous")
gg <- Hmisc::cut2(geneCount, g=ng)
plot(rd, pch=16, cex=1/2, col=pal[gg], main=main, ...)
}
# progenitors
plotByGene(rd[, 1:2], countsFQ["Dpp4",], ng=2, main="Dpp4")
plotByGene(rd[, 1:2], countsFQ["Wnt2",], ng=4, main="Wnt2")
# group 2 cells
plotByGene(rd[, 1:2], countsFQ["Icam1",], ng=4, main="Icam1")
plotByGene(rd[, 1:2], countsFQ["Dlk1",], ng=4, main="Pref1") #Pref1 is also called Dlk1
# group 3 cells
plotByGene(rd[, 1:2], countsFQ["Clec11a",], ng=4, main="Clec11a")
# group 4 cells
plotByGene(rd[, 1:2], countsFull["Wnt6",], ng=2, main="Wnt6")
# cluster
set.seed(9)
cl <- kmeans(rd, centers = 10)
color = grDevices::colors()[grep('gr(a|e)y', grDevices::colors(), invert = T)]
pal <- sample(color, 12)
plot(rd[,1:2], pch=16, cex=1/2, col=pal[cl$cluster])
legend("bottomleft", legend=1:10, pch=16, col=pal, bty='n')
# cluster 3 and 7 are progenitor cells
# cluster 1, 5 and 10 are group 2 cells
# cluster 2 and 6 are group 3 cells
keepCells <- which(cl$cluster %in% c(1,2,3,5,7,10))
countsFQ2 <- countsFQ[,keepCells]
rd2 <- prcomp(t(log1p(countsFQ2)), scale. = FALSE)
### umap
library(umap)
rdUmap <- umap(rd2$x[,1:20])
plot(rdUmap$layout, pch=16, cex=1/3)
rafalib::mypar(mfrow=c(3,2))
# progenitors
plotByGene(rdUmap$layout, countsFQ2["Dpp4",], ng=2, main="Prog: Dpp4")
plotByGene(rdUmap$layout, countsFQ2["Wnt2",], ng=4, main="Prog: Wnt2")
# group 2 cells
plotByGene(rdUmap$layout, countsFQ2["Icam1",], ng=4, main="G2: Icam1")
plotByGene(rdUmap$layout, countsFQ2["Dlk1",], ng=4, main="G2: Pref1") #Pref1 is also called Dlk1
# group 3 cells
plotByGene(rdUmap$layout, countsFQ2["Clec11a",], ng=4, main="G3: Clec11a")
# group 4 cells
plotByGene(rdUmap$layout, countsFull["Wnt6",keepCells], ng=2, main="G4: Wnt6")
# cluster
set.seed(12)
pal <- wesanderson::wes_palette("Darjeeling1", n=6, type="continuous")
cl <- kmeans(rdUmap$layout, centers = 6)
plot(rdUmap$layout, pch=16, cex=1/2, col=pal[cl$cluster])
legend("topleft", legend=1:6, pch=16, col=pal[1:6], bty='n')
# slingshot
library(slingshot)
## Loading required package: princurve
lin <- getLineages(rdUmap$layout, clusterLabels=cl$cluster,
start.clus=1, end.clus=c(5,2))
## Using full covariance matrix
plot(rdUmap$layout, pch=16, cex=1/2, col=pal[cl$cluster])
lines(lin, lwd=2)
crv <- getCurves(lin)
## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values
## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values
## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values
## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values
## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values
## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values
## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values
## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values
## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values
## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values
## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values
## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values
## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values
## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values
## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values
## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values
## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values
## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values
## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values
## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values
## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values
## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values
## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values
## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values
## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values
## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values
## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values
## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values
## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values
## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values
## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values
## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values
## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values
## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values
## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values
## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values
## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values
## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values
## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values
## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values
## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values
## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values
## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values
## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values
## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values
## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values
## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values
## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values
## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values
## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values
## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values
## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values
## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values
## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values
## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values
## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values
## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values
## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values
## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values
## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values
## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values
## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values
## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values
## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values
## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values
## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values
## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values
## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values
## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values
## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values
## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values
## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values
## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values
## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values
## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values
## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values
## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values
## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values
## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values
## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values
## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values
## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values
## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values
## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values
## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values
## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values
## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values
## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values
## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values
## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values
## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values
## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values
## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values
## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values
## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values
## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values
## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values
## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values
## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values
## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values
## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values
## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values
## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values
## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values
## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values
## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values
## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values
## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values
## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values
## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values
## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values
## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values
## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values
## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values
## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values
## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values
## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values
## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values
## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values
## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values
plot(rdUmap$layout, pch=16, cex=1/2, col=pal[cl$cluster],
xlab="UMAP1", ylab="UMAP2", bty='l')
lines(crv, lwd=2, col="black")
saveRDS(crv, file="crvUmap.rds")
8 knots seem appropriate.
library(tradeSeq)
## tradeSeq has been updated to accommodate singleCellExperiment objects as output, making it much more memory efficient. Please check the news file and the updated vignette for details.
set.seed(3)
aicMat1 <- evaluateK(countsFQ2, k=3:10, nGenes=200, sds=crv)
set.seed(4)
aicMat2 <- evaluateK(countsFQ2, k=3:10, nGenes=200, sds=crv)
sce <- fitGAM(countsFQ2, sds=crv, nknots=8)
saveRDS(sce, file="~/data/sceAdipose.rda")
seRes <- startVsEndTest(sce)
ose <- order(seRes$waldStat, decreasing=TRUE)
which(rownames(seRes)[ose] == "Wnt2") # ranked 263
## [1] 232
which(rownames(seRes)[ose] == "Dpp4") # ranked 450
## [1] 449
head(seRes[ose,], n=16)
## waldStat df pvalue
## Igfbp7 2360.6550 2 0
## Mfap5 1718.0559 2 0
## Clec3b 1563.6112 2 0
## Pi16 1557.9152 2 0
## Apoe 1324.7050 2 0
## Akr1c18 1296.7321 2 0
## Igfbp6 1196.7044 2 0
## Fn1 1181.2875 2 0
## Anxa1 1118.8296 2 0
## Fbn1 1107.9533 2 0
## Bgn 1098.3094 2 0
## Ly6a 1023.9972 2 0
## Gap43 985.6586 2 0
## Fstl1 941.2143 2 0
## Ly6c1 875.0512 2 0
## Pcolce2 777.1175 2 0
rafalib::mypar(mfrow=c(3,2))
sapply(1:6, function(ii){
gene <- rownames(seRes)[ose][ii]
print(plotByGene(rdUmap$layout, countsFQ2[gene,], ng=4, main=gene))
})
## NULL
## NULL
## NULL
## NULL
## NULL
## NULL
## [[1]]
## NULL
##
## [[2]]
## NULL
##
## [[3]]
## NULL
##
## [[4]]
## NULL
##
## [[5]]
## NULL
##
## [[6]]
## NULL
sapply(7:12, function(ii){
gene <- rownames(seRes)[ose][ii]
print(plotByGene(rdUmap$layout, countsFQ2[gene,], ng=4, main=gene))
})
## NULL
## NULL
## NULL
## NULL
## NULL
## NULL
## [[1]]
## NULL
##
## [[2]]
## NULL
##
## [[3]]
## NULL
##
## [[4]]
## NULL
##
## [[5]]
## NULL
##
## [[6]]
## NULL
#figure for paper
rafalib::mypar(mfrow=c(3,2))
genes <- c("Dpp4", "Wnt2",
"Pi16", "Akr1c18",
"Fn1", "Fbn1")
sapply(genes, function(gene){
print(plotByGene(rdUmap$layout, countsFQ2[gene,], ng=4, main=gene, bty='n',
xlab="UMAP1", ylab="UMAP2"))
})
## NULL
## NULL
## NULL
## NULL
## NULL
## NULL
## $Dpp4
## NULL
##
## $Wnt2
## NULL
##
## $Pi16
## NULL
##
## $Akr1c18
## NULL
##
## $Fn1
## NULL
##
## $Fbn1
## NULL
Note that this is only restricted to the end points and is therefore limited.
deRes <- diffEndTest(sce)
ode <- order(deRes$waldStat, decreasing=TRUE)
head(seRes[ode,], n=10)
## waldStat df pvalue
## Meox2 388.2344 2 0
## Igfbp7 2360.6550 2 0
## Cst3 525.6228 2 0
## Nr4a1 203.5429 2 0
## Igfbp6 1196.7044 2 0
## Tmsb4x 712.0800 2 0
## Col14a1 702.0305 2 0
## Cpe 247.7400 2 0
## Cd81 244.5232 2 0
## Bgn 1098.3094 2 0
rafalib::mypar(mfrow=c(3,2))
sapply(1:6, function(ii){
gene <- rownames(deRes)[ode][ii]
print(plotByGene(rdUmap$layout, countsFQ2[gene,], ng=4, main=gene))
})
## NULL
## NULL
## NULL
## NULL
## NULL
## NULL
## [[1]]
## NULL
##
## [[2]]
## NULL
##
## [[3]]
## NULL
##
## [[4]]
## NULL
##
## [[5]]
## NULL
##
## [[6]]
## NULL
plotGeneCount(crv, countsFQ2, models=sce, clusters = cl$cluster)
resedt <- earlyDETest(sce, knots=c(3,7))
oedt <- order(resedt$waldStat, decreasing=TRUE)
head(resedt[oedt,], n=10)
## waldStat df pvalue
## Dlk1 3521.369 7 0
## H19 3028.316 7 0
## Itm2a 2775.809 7 0
## Mgp 2524.158 7 0
## Rarres2 2436.242 7 0
## Meox2 2188.965 7 0
## Postn 2085.855 7 0
## Fabp4 2062.955 7 0
## Cst3 1934.916 7 0
## Igfbp6 1766.272 7 0
rafalib::mypar(mfrow=c(3,2))
sapply(1:6, function(ii){
gene <- rownames(resedt)[oedt][ii]
print(plotByGene(rdUmap$layout, countsFQ2[gene,], ng=4, main=gene))
})
## NULL
## NULL
## NULL
## NULL
## NULL
## NULL
## [[1]]
## NULL
##
## [[2]]
## NULL
##
## [[3]]
## NULL
##
## [[4]]
## NULL
##
## [[5]]
## NULL
##
## [[6]]
## NULL
#figure for paper
rafalib::mypar(mfrow=c(1,2))
genes <- c("Mgp", "Meox2")
sapply(genes, function(gene){
print(plotByGene(rdUmap$layout, countsFQ2[gene,], ng=4, main=gene, bty='n',
xlab="UMAP1", ylab="UMAP2"))
})
## NULL
## NULL
## $Mgp
## NULL
##
## $Meox2
## NULL
#figure for paper
rafalib::mypar(mfrow=c(1,2))
genes <- c("H19", "Col14a1")
sapply(genes, function(gene){
print(plotByGene(rdUmap$layout, countsFQ2[gene,], ng=4, main=gene, bty='n',
xlab="UMAP1", ylab="UMAP2"))
})
## NULL
## NULL
## $H19
## NULL
##
## $Col14a1
## NULL